home *** CD-ROM | disk | FTP | other *** search
/ Practical Algorithms for Image Analysis / Practical Algorithms for Image Analysis.iso / CH_4.3 / XCP / BDY.C next >
Encoding:
C/C++ Source or Header  |  1999-09-11  |  15.1 KB  |  599 lines
  1. /* 
  2.  * bdy.c
  3.  * 
  4.  * Practical Algorithms for Image Analysis
  5.  * 
  6.  * Copyright (c) 1997, 1998, 1999 MLMSoftwareGroup, LLC
  7.  */
  8.  
  9. /*
  10.  * BDY.C
  11.  *
  12.  * functions implementating C.T.Zahn's algorithm for boundary
  13.  * detection and representation
  14.  *
  15.  */
  16. #include <stdio.h>
  17. #include <malloc.h>
  18. #include <stdlib.h>
  19. #include <math.h>
  20. #include "ip.h"
  21. #include "bdy.h"
  22.  
  23. #define    WINDOW_ROWS    3
  24. #define NDIRECTIONS    8
  25.  
  26. #define    THRESH        1
  27. #define MIN_POLY_SIZE    1
  28. #define    ZERO        0
  29. #define ROOT2        1.414213562
  30.  
  31. #define    ON        1
  32. #define    OFF        0
  33. #undef    SHOW_CURV_PT
  34.  
  35. #define    RESET        ON
  36.  
  37. #define IN_LIST        1
  38. #define    OUT_LIST    0
  39. #define SHOW_CURV_PT
  40.  
  41.  
  42. /* global variables */
  43. extern struct curv_points *curv_head_in[NDIRECTIONS];
  44. extern struct curv_points *curv_head_out[NDIRECTIONS];
  45. extern struct curv_points *curv_tail_in[NDIRECTIONS];
  46. extern struct curv_points *curv_tail_out[NDIRECTIONS];
  47.  
  48. extern struct polygon *poly_head_ptr;
  49. extern struct polygon *poly_tail_ptr;
  50.  
  51. extern float *delta_phik, *delta_lk;
  52. extern long ncp;
  53.  
  54. unsigned int d1, d2, d3, d4, d5, d6, b1, b2, b3, b4, b5, b6;
  55. unsigned nbytes = 1;
  56.  
  57.  
  58. /*
  59.  * generate curvature points using Zahn algorithm
  60.  * (ref: C. T. Zahn, SLAC report No. 70, Dec. 1966)
  61.  */
  62. void
  63. curvature_points (unsigned char *window[], int ir, int jmax, int imax)
  64. {
  65.   int jc;
  66.  
  67.   for (jc = 0; jc < jmax - 2; jc++) {
  68.     d1 = *(*(window + (ir % WINDOW_ROWS)) + jc);
  69.     d2 = *(*(window + (ir % WINDOW_ROWS)) + jc + 1);
  70.     d3 = *(*(window + (ir % WINDOW_ROWS)) + jc + 2);
  71.     d4 = *(*(window + ((ir + 1) % WINDOW_ROWS)) + jc);
  72.     d5 = *(*(window + ((ir + 1) % WINDOW_ROWS)) + jc + 1);
  73.     d6 = *(*(window + ((ir + 1) % WINDOW_ROWS)) + jc + 2);
  74.     horizontal (ir, jc);
  75.   }
  76.   for (jc = 0; jc < jmax - 1; jc++) {
  77.     b1 = *(*(window + (ir % WINDOW_ROWS)) + jc + 1);
  78.     b2 = *(*(window + ((ir + 1) % WINDOW_ROWS)) + jc + 1);
  79.     b3 = *(*(window + ((ir + 2) % WINDOW_ROWS)) + jc + 1);
  80.     b4 = *(*(window + (ir % WINDOW_ROWS)) + jc);
  81.     b5 = *(*(window + ((ir + 1) % WINDOW_ROWS)) + jc);
  82.     b6 = *(*(window + ((ir + 2) % WINDOW_ROWS)) + jc);
  83.     vertical (ir, jc);
  84.   }
  85.   if (ir == imax - WINDOW_ROWS) {
  86.     ir++;
  87.     for (jc = 0; jc < jmax - 2; jc++) {
  88.       d1 = *(*(window + (ir % WINDOW_ROWS)) + jc);
  89.       d2 = *(*(window + (ir % WINDOW_ROWS)) + jc + 1);
  90.       d3 = *(*(window + (ir % WINDOW_ROWS)) + jc + 2);
  91.       d4 = *(*(window + ((ir + 1) % WINDOW_ROWS)) + jc);
  92.       d5 = *(*(window + ((ir + 1) % WINDOW_ROWS)) + jc + 1);
  93.       d6 = *(*(window + ((ir + 1) % WINDOW_ROWS)) + jc + 2);
  94.       horizontal (ir, jc);
  95.     }
  96.   }
  97. }
  98.  
  99.  
  100. /*
  101.  * check 3x2 horizontal array (in a binary image) and determine 
  102.  * direction codes for incoming and outgoing contour segment
  103.  */
  104. int
  105. horizontal (i, j)
  106.      int i, j;
  107. {
  108.   int in, out;
  109.   int err_flag = -1;
  110.  
  111.   if (d2 >= THRESH) {
  112.     if (d5 == ZERO) {
  113.       if (d6 >= THRESH)
  114.         in = 3;
  115.       else if (d3 >= THRESH)
  116.         in = 4;
  117.       else
  118.         in = 5;
  119.  
  120.       if (d4 >= THRESH)
  121.         out = 5;
  122.       else if (d1 >= THRESH)
  123.         out = 4;
  124.       else
  125.         out = 3;
  126.     }
  127.     else
  128.       return (err_flag);
  129.   }
  130.   else {
  131.     if (d5 >= THRESH) {
  132.       if (d1 >= THRESH)
  133.         in = 7;
  134.       else if (d4 >= THRESH)
  135.         in = 0;
  136.       else
  137.         in = 1;
  138.  
  139.       if (d3 >= THRESH)
  140.         out = 1;
  141.       else if (d6 >= THRESH)
  142.         out = 0;
  143.       else
  144.         out = 7;
  145.     }
  146.     else
  147.       return (err_flag);
  148.   }
  149.  
  150. /*
  151.  * curvature point found! 
  152.  * allocate memory and store it in linked list curv_points structure 
  153.  */
  154.   if (in != out) {
  155.  
  156. #ifdef SHOW_CURV_PT
  157.     printf ("   horizontal:  ir= %d, jc= %d, in= %d, out= %d\n",
  158.             2 * (i + 1), 2 * (j + 1) + 1, in, out);
  159. #endif
  160.     if (curv_head_in[in] == NULL) {
  161.       curv_head_in[in] = (struct curv_points *) calloc (nbytes, sizeof (struct curv_points));
  162.       curv_head_in[in]->prev = NULL;
  163.       curv_tail_in[in] = curv_head_in[in];
  164.     }
  165.     else {
  166.       curv_tail_in[in]->next = (struct curv_points *) calloc (nbytes, sizeof (struct curv_points));
  167.       curv_tail_in[in]->next->prev = curv_tail_in[in];
  168.       curv_tail_in[in] = curv_tail_in[in]->next;
  169.     }
  170.     if (curv_head_out[out] == NULL) {
  171.       curv_head_out[out] = (struct curv_points *) calloc (nbytes, sizeof (struct curv_points));
  172.       curv_head_out[out]->prev = NULL;
  173.       curv_tail_out[out] = curv_head_out[out];
  174.     }
  175.     else {
  176.       curv_tail_out[out]->next = (struct curv_points *) calloc (nbytes, sizeof (struct curv_points));
  177.       curv_tail_out[out]->next->prev = curv_tail_out[out];
  178.       curv_tail_out[out] = curv_tail_out[out]->next;
  179.     }
  180.     ncp++;
  181.  
  182.     curv_tail_in[in]->same_point = curv_tail_out[out];
  183.     curv_tail_in[in]->next = NULL;
  184.     curv_tail_out[out]->next = NULL;
  185.     curv_tail_in[in]->xn = curv_tail_out[out]->xn = 2 * (j + 1) + 1;
  186.     curv_tail_in[in]->yn = curv_tail_out[out]->yn = 2 * (i + 1);
  187.     curv_tail_in[in]->edge_in = curv_tail_out[out]->edge_in = in;
  188.     curv_tail_in[in]->edge_out = curv_tail_out[out]->edge_out = out;
  189.   }
  190.   return (1);
  191. }
  192.  
  193.  
  194. /*
  195.  * check 2x3 vertical array (in a binary image) and determine 
  196.  * direction codes for incoming and outgoing contour segment
  197.  */
  198. int
  199. vertical (i, j)
  200.      int i, j;
  201. {
  202.   int in, out;
  203.   int err_flag = -1;
  204.  
  205.   if (b2 >= THRESH) {
  206.     if (b5 == ZERO) {
  207.       if (b6 >= THRESH)
  208.         in = 1;
  209.       else if (b3 >= THRESH)
  210.         in = 2;
  211.       else
  212.         in = 3;
  213.  
  214.       if (b4 >= THRESH)
  215.         out = 3;
  216.       else if (b1 >= THRESH)
  217.         out = 2;
  218.       else
  219.         out = 1;
  220.     }
  221.     else
  222.       return (err_flag);
  223.   }
  224.   else {
  225.     if (b5 >= THRESH) {
  226.       if (b1 >= THRESH)
  227.         in = 5;
  228.       else if (b4 >= THRESH)
  229.         in = 6;
  230.       else
  231.         in = 7;
  232.  
  233.       if (b3 >= THRESH)
  234.         out = 7;
  235.       else if (b6 >= THRESH)
  236.         out = 6;
  237.       else
  238.         out = 5;
  239.     }
  240.     else
  241.       return (err_flag);
  242.   }
  243.  
  244. /*
  245.  * curvature point found! 
  246.  * allocate memory and store it in linked list curv_points structure 
  247.  */
  248.   if (in != out) {
  249.  
  250. #ifdef SHOW_CURV_PT
  251.     printf ("  vertical:  ir= %d, jc= %d, in= %d, out= %d\n",
  252.             2 * (i + 1) + 1, 2 * (j + 1), in, out);
  253. #endif
  254.  
  255.     if (curv_head_in[in] == NULL) {
  256.       curv_head_in[in] = (struct curv_points *) calloc (nbytes, sizeof (struct curv_points));
  257.       curv_head_in[in]->prev = NULL;
  258.       curv_tail_in[in] = curv_head_in[in];
  259.     }
  260.     else {
  261.       curv_tail_in[in]->next = (struct curv_points *) calloc (nbytes, sizeof (struct curv_points));
  262.       curv_tail_in[in]->next->prev = curv_tail_in[in];
  263.       curv_tail_in[in] = curv_tail_in[in]->next;
  264.     }
  265.     if (curv_head_out[out] == NULL) {
  266.       curv_head_out[out] = (struct curv_points *) calloc (nbytes, sizeof (struct curv_points));
  267.       curv_head_out[out]->prev = NULL;
  268.       curv_tail_out[out] = curv_head_out[out];
  269.     }
  270.     else {
  271.       curv_tail_out[out]->next = (struct curv_points *) calloc (nbytes, sizeof (struct curv_points));
  272.       curv_tail_out[out]->next->prev = curv_tail_out[out];
  273.       curv_tail_out[out] = curv_tail_out[out]->next;
  274.     }
  275.     ncp++;
  276.  
  277.     curv_tail_in[in]->same_point = curv_tail_out[out];
  278.     curv_tail_in[in]->next = NULL;
  279.     curv_tail_out[out]->next = NULL;
  280.     curv_tail_in[in]->xn = curv_tail_out[out]->xn = 2 * (j + 1);
  281.     curv_tail_in[in]->yn = curv_tail_out[out]->yn = 2 * (i + 1) + 1;
  282.     curv_tail_in[in]->edge_in = curv_tail_out[out]->edge_in = in;
  283.     curv_tail_in[in]->edge_out = curv_tail_out[out]->edge_out = out;
  284.   }
  285.   return (1);
  286. }
  287.  
  288.  
  289. /*
  290.  * link up all curvature points into closed polygons
  291.  */
  292. void
  293. linkage (Image * imgIn, Image * imgOut, int value)
  294. {
  295.   struct polygon *dummy_poly_ptr;
  296.   struct curv_points *first_point;
  297.   float *d_phi_head, *d_l_head;
  298.   int m, c;
  299.   char inbuf[IN_BUF_LEN];
  300.  
  301.   do {
  302.     d_phi_head = delta_phik;
  303.     d_l_head = delta_lk;
  304.     for (m = 0; m < NDIRECTIONS; m++)
  305.       if ((first_point = curv_head_out[m]) != NULL) {
  306.         poly (first_point);
  307.         break;
  308.       }
  309.   } while (first_point != NULL);
  310.  
  311.   if (ncp != 0)
  312.     printf ("...something wrong, didn't find all polygons!\n");
  313.   if (poly_head_ptr == NULL) {
  314.     printf ("...no objects found to analyze!\n");
  315.     return;
  316.   }
  317.  
  318. /*
  319.  * polygon analysis
  320.  */
  321.   printf ("...generate curvature points for polygons?(y/n) -- ");
  322.   while ((c = readlin (inbuf)) == 'y') {
  323.     dummy_poly_ptr = select_poly (poly_head_ptr, imgOut, value);
  324.     printf ("\n...generate curvature points for another polygon?(y/n)");
  325.   }
  326. }
  327.  
  328. /*
  329.  * create a closed polygon given a starting curvature point 
  330.  */
  331. void
  332. poly (first_point)
  333.      struct curv_points *first_point;
  334. {
  335.   struct curv_points *next_point, *prev_point, *match_in_list ();
  336.  
  337.   if (poly_head_ptr == NULL) {
  338.     poly_head_ptr = (struct polygon *) calloc (nbytes, sizeof (struct polygon));
  339.     poly_tail_ptr = poly_head_ptr;
  340.     poly_head_ptr->prev_poly = NULL;
  341.   }
  342.   else {
  343.     poly_tail_ptr->next_poly = (struct polygon *) calloc (nbytes, sizeof (struct polygon));
  344.     poly_tail_ptr->next_poly->prev_poly = poly_tail_ptr;
  345.     poly_tail_ptr = poly_tail_ptr->next_poly;
  346.   }
  347.   poly_tail_ptr->next_poly = NULL;
  348.  
  349.   prev_point = first_point;
  350.   next_point = match_in_list (prev_point);
  351.   poly_tail_ptr->d_phi_ptr = delta_phik;
  352.   poly_tail_ptr->d_l_ptr = delta_lk;
  353.   poly_tail_ptr->poly_points = 0;
  354.   poly_tail_ptr->total_delta_phi = (float) 0;
  355.   poly_tail_ptr->first_x = first_point->xn;
  356.   poly_tail_ptr->first_y = first_point->yn;
  357.   poly_tail_ptr->first_edge_out = first_point->edge_out;
  358.   create_edge (prev_point, next_point);
  359.   delete_list (prev_point, OUT_LIST);
  360.   prev_point = next_point->same_point;
  361.   delete_list (next_point, IN_LIST);
  362.   ncp--;
  363.   poly_tail_ptr->poly_points++;
  364.   poly_tail_ptr->total_delta_phi = poly_tail_ptr->total_delta_phi + *delta_phik;
  365.   delta_lk++;
  366.   delta_phik++;
  367.   do {
  368.     next_point = match_in_list (prev_point);
  369.     create_edge (prev_point, next_point);
  370.     delete_list (prev_point, OUT_LIST);
  371.     prev_point = next_point->same_point;
  372.     delete_list (next_point, IN_LIST);
  373.     ncp--;
  374.     poly_tail_ptr->poly_points++;
  375.     poly_tail_ptr->total_delta_phi = poly_tail_ptr->total_delta_phi
  376.       + *delta_phik;
  377.     delta_lk++;
  378.     delta_phik++;
  379.  
  380.     if (ncp < 0) {
  381.       printf ("Something wrong, points<0\n");
  382.       exit (1);
  383.     }
  384.   }
  385.   while (prev_point != first_point);
  386.  
  387.   if (poly_tail_ptr->poly_points <= MIN_POLY_SIZE) {
  388.     if (poly_tail_ptr->prev_poly == NULL)
  389.       poly_head_ptr = NULL;
  390.     else {
  391.       poly_tail_ptr = poly_tail_ptr->prev_poly;
  392.       free (poly_tail_ptr->next_poly);
  393.       poly_tail_ptr->next_poly = NULL;
  394.     }
  395.   }
  396. }
  397.  
  398. /*
  399.  * given a point in the out-list, find the next point
  400.  * in the in-list with edge-out(out-list) = edge-in(in-list)
  401.  */
  402. struct curv_points *
  403. match_in_list (current_point)
  404.      struct curv_points *current_point;
  405. {
  406.   struct curv_points *tail_point, *head_point;
  407.   int direction;
  408.   int xcurr, ycurr, xnext, ynext;
  409.  
  410.   xcurr = current_point->xn;
  411.   ycurr = current_point->yn;
  412.   printf ("match_in_list: current point = (%d, %d)\n", xcurr, ycurr);
  413.   direction = current_point->edge_out;
  414.   tail_point = curv_tail_in[direction];
  415.   head_point = curv_head_in[direction];
  416.   while ((head_point != NULL) || (tail_point != NULL)) {
  417.     switch (direction) {
  418.     case 0:
  419.       ynext = head_point->yn;
  420.       xnext = head_point->xn;
  421.       if ((ycurr == ynext) && (xnext > xcurr))
  422.         return (head_point);
  423.       break;
  424.     case 4:
  425.       ynext = tail_point->yn;
  426.       xnext = tail_point->xn;
  427.       if ((ycurr == ynext) && (xnext < xcurr))
  428.         return (tail_point);
  429.       break;
  430.     case 2:
  431.       ynext = tail_point->yn;
  432.       xnext = tail_point->xn;
  433.       if ((xcurr == xnext) && (ynext < ycurr))
  434.         return (tail_point);
  435.       break;
  436.     case 6:
  437.       ynext = head_point->yn;
  438.       xnext = head_point->xn;
  439.       if ((xcurr == xnext) && (ynext > ycurr))
  440.         return (head_point);
  441.       break;
  442.     case 1:
  443.       ynext = tail_point->yn;
  444.       xnext = tail_point->xn;
  445.       if (((xcurr + ycurr) == (xnext + ynext)) &&
  446.           (xnext > xcurr) && (ynext < ycurr))
  447.         return (tail_point);
  448.       break;
  449.     case 5:
  450.       ynext = head_point->yn;
  451.       xnext = head_point->xn;
  452.       if (((xcurr + ycurr) == (xnext + ynext)) &&
  453.           (xnext < xcurr) && (ynext > ycurr))
  454.         return (head_point);
  455.       break;
  456.     case 3:
  457.       ynext = tail_point->yn;
  458.       xnext = tail_point->xn;
  459.       if (((xcurr - ycurr) == (xnext - ynext)) &&
  460.           (xnext < xcurr) && (ynext < ycurr))
  461.         return (tail_point);
  462.       break;
  463.     case 7:
  464.       ynext = head_point->yn;
  465.       xnext = head_point->xn;
  466.       if (((xcurr - ycurr) == (xnext - ynext)) &&
  467.           (xnext > xcurr) && (ynext > ycurr))
  468.         return (head_point);
  469.       break;
  470.     }
  471.     tail_point = tail_point->prev;
  472.     head_point = head_point->next;
  473.   }
  474.   printf ("...no match found in outlist[%d]\n", direction);
  475.   exit (1);
  476. }
  477.  
  478. /*
  479.  * delete curvature point from either in-list or out-list
  480.  */
  481. void
  482. delete_list (current_point, list)
  483.      struct curv_points *current_point;
  484.      int list;
  485. {
  486.   int direction;
  487.  
  488.   if (current_point == NULL)
  489.     printf ("...current point can't be NULL!\n");
  490.  
  491.   if ((current_point->prev == NULL) && (current_point->next == NULL)) {
  492.     if (list == IN_LIST) {
  493.       direction = current_point->edge_in;
  494.       curv_head_in[direction] = NULL;
  495.       curv_tail_in[direction] = NULL;
  496.     }
  497.     else {
  498.       direction = current_point->edge_out;
  499.       curv_head_out[direction] = NULL;
  500.       curv_tail_out[direction] = NULL;
  501.     }
  502.   }
  503.   else if (current_point->prev == NULL) {
  504.     if (list == IN_LIST) {
  505.       direction = current_point->edge_in;
  506.       curv_head_in[direction] = current_point->next;
  507.     }
  508.     else {
  509.       direction = current_point->edge_out;
  510.       curv_head_out[direction] = current_point->next;
  511.     }
  512.     current_point->next->prev = NULL;
  513.   }
  514.   else if (current_point->next == NULL) {
  515.     if (list == IN_LIST) {
  516.       direction = current_point->edge_in;
  517.       curv_tail_in[direction] = current_point->prev;
  518.     }
  519.     else {
  520.       direction = current_point->edge_out;
  521.       curv_tail_out[direction] = current_point->prev;
  522.     }
  523.     current_point->prev->next = NULL;
  524.   }
  525.   else {
  526.     current_point->prev->next = current_point->next;
  527.     current_point->next->prev = current_point->prev;
  528.   }
  529.   free (current_point);
  530. }
  531.  
  532.  
  533. void
  534. print_curv ()
  535. {
  536.   int x, y;
  537.   struct curv_points *temp_head;
  538.  
  539.   for (x = 0; x < NDIRECTIONS; x++) {
  540.     y = 0;
  541.     temp_head = curv_head_in[x];
  542.     while (temp_head != NULL) {
  543.       y++;
  544.       temp_head = temp_head->next;
  545.     }
  546.     if (y != 0)
  547.       printf ("inlist[%d] = %d\n", x, y);
  548.   }
  549.  
  550.   for (x = 0; x < NDIRECTIONS; x++) {
  551.     y = 0;
  552.     temp_head = curv_head_out[x];
  553.     while (temp_head != NULL) {
  554.       y++;
  555.       temp_head = temp_head->next;
  556.     }
  557.     if (y != 0)
  558.       printf ("outlist[%d] = %d\n", x, y);
  559.   }
  560. }
  561.  
  562.  
  563. /*
  564.  * create an edge from old-point to new-point
  565.  */
  566. void
  567. create_edge (old_point, new_point)
  568.      struct curv_points *old_point, *new_point;
  569. {
  570.   int phi;
  571.  
  572.   switch (new_point->edge_in) {
  573.   case 0:
  574.   case 4:
  575.     *delta_lk = (float) abs (new_point->xn - old_point->xn);
  576.     phi = new_point->edge_out - new_point->edge_in;
  577.     if (phi == 7)
  578.       phi = -1;
  579.     *delta_phik = (float) phi;
  580.     break;
  581.   case 2:
  582.   case 6:
  583.     *delta_lk = (float) abs (new_point->yn - old_point->yn);
  584.     *delta_phik = (float) (new_point->edge_out - new_point->edge_in);
  585.     break;
  586.   default:
  587.     *delta_lk = (float) abs (new_point->xn - old_point->xn) * (float) ROOT2;
  588.     phi = new_point->edge_out - new_point->edge_in;
  589.     if (phi == -7)
  590.       phi = 1;
  591.     else if (phi == -6)
  592.       phi = 2;
  593.     else if (phi == 6)
  594.       phi = -2;
  595.     *delta_phik = (float) phi;
  596.     break;
  597.   }
  598. }
  599.